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ABSTRACT 



We present an efficient solution to the following problem, of 
relevance in a numerical optimization scheme: calculation of 
integrals of the type 

01 02 dxdy 
Tn{/>0} 

for quadratic polynomials /, 01,02 on a plane triangle T. 
The naive approach would involve consideration of the many 
possible shapes of T fl {/ > 0} (possibly after a convenient 
transformation) and parameterizing its border, in order to 
integrate the variables separately. Our solution involves par- 
titioning the triangle into smaller triangles on which integra- 
tion is much simpler. 

Categories and Subject Descriptors 

G.1.8 [Numerical Analysis]: Partial Differential Equa- 
tions — finite element methods; 1.1.2 [Symbolic and alge- 
braic manipulation]: Algorithms — algebraic algorithms 



Keywords 

symbolic integration, triangular subdivision, optimal con- 
trol, variational discretization, quadratic shape functions 



1. INTRODUCTION 

This article presents a symbolic solution to a problem of 
relevance in a numerical optimization scheme: the numerical 
solution of optimal control problems with partial differential 
equations as constraint requires to discretize the problem, 
i.e. to solve finite-dimensional approximations, see e.g. [^. 
When applying the variational discretization concept [3] , the 
following problem arises: integrals of the type 



//, 



Tn{/>0} 



g{x,y)dxdy, 



9 = (l>i 



(1) 
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for quadratic polynomials /, 0i , 02 on a plane triangle T have 
to be evaluated accurately. Up to now the variational dis- 
cretization method was used only for degree 1, i.e. where 
the function / defining the integral region in ([1} is a polyno- 
mial of degree 1. Using polynomials with higher order gives 
better approximation results, see Theorem 15.5! below. 

The naive approach to compute the integral ^ would 
involve consideration of the many possible shapes of rn{/ > 
0} and parameterizing its border, in order to integrate the 
variables separately. This suffers from some computational 
difficulties as we show below. 

Example 1.1. Suppose that one side of the triangle lies 
on a horizontal line. Consider the situation where the region 
of integration is the part of the interior of an ellipse in the 
triangle, as in the figure. 




Xl X2'X3 -x-i 



We see that we could calculate the integral as the sum of 
five integrals on domains perpendicular to the x axis. For 
example the first one is 



-(^) 



9{x,y)dy dx 



where y — Iac{^) is the equation of the line AC and y = 
c_ (a;) is the equation of the lower part of the ellipse. Note 
that we need to parameterize the ellipse (this will involve at 
best a square root or trigonometric functions) and calculate 
the x-coordmates of the relevant points, which also involve 
square roots. The value of the inner integral will be given by 
a formula of which an antiderivative must be computed then. 
The resulting formula is far from simple. 

An alternative would be to apply an affine transformation 
so that the ellipse becomes a circle centered at the origin, 
followed by a change to polar coordinates. This does not 
make the integral significantly easier to compute. 

And of course, here we use our knowledge of the relative 
position of the ellipse and the triangle, as in the figure; the 



possible relative positions of a conic and a triangle are many, 
and to discern them is not trivial. 

In contrast, in the following particular case we obtain a 
simple formula. 

Example 1.2. Let A = (0, 0), B = (1, 0), C = (1, 1), and 
assume that f > on the triangle T — ABC. If g{x,y) — 
y^ bijX^y-' then the integral becomes 



IL 



g{x,y)dxdy= j i g{x,y)dy\dx^ 

= E r. , J:\ . , o^ (2) 



i+j<4 



iJ + l){i + j + 2) 



For a general triangle T, one applies an affine transfor- 
mation which brings the vertices to the above points. The 
resulting integrand is a polynomial of the same degree, and 
only a constant factor is introduced by the substitution for- 
mula. 

Our solution involves partitioning T into smaller triangles 
on which integration is much simpler. The result is a decision 
tree and several relatively simple explicit formulas, which 
form Algorithm 14.11 The particular nature of g, beyond it 
being a polynomial, will be immaterial. Our method could 
in principle be adapted for larger values of deg /, although 
it may become too complicated for practical uses even for 
degree 3. Besides, only the quadratic case is relevant for 
the context in which this problem arose. It is important to 
point out that an implementation for floating point arith- 
metic would need a more detailed treatment, see the end of 
Section [S] 

We describe the subdivision method in Section (2] the in- 
tegration in the base cases in Section O the complete algo- 
rithm in Section |4l and a description of the application to 
optimization in Section [S] which includes some comments on 
the practical implementation of our algorithm. 

2. TRIANGULAR SUBDIVISION 

Our idea is to reduce the number of intersections between 
the curve / = and the sides of the triangle, by cutting 
the triangle into pieces until we reach some base cases that 
we establish below. For those cases the integration will be 
much simpler than in Example 11.11 We leave for later the 
case where the conic / = is degenerate (two lines, either 
intersecting, parallel, or coincident; one point; the empty 
set). Note that the type of a conic can be determined quickly 
by inspection of the equation. 

First we introduce some nomenclature. 

Definition 2.1. Fix a nonsmgular conic. A segment is 
called free (with respect to the conic) if it does not intersect 
it except possibly at the vertices of the segment. 

Remark 2.2. 

1. Any line or segment intersects any conic at most at 
two points. 

2. A segment joining two points of the conic is always 
free. 



Proof. Part 1 is a simple case of the weak Bezout's the- 
orem [5]. Part 2 is a clear consequence of part 1. D 

The calculation of the intersections of a segment with a 
given conic (and thus the determination of the freedom of 
the segment) is straightforward, and fast in practice. 

The next definition encapsulates the base cases of our sub- 
division method. 

Definition 2.3. A triangle is called free if either all its 
sides are free, or the intersection of its border with the conic 
is just one non-vertex point. 

Thus there are five types of free triangles: those with all 
sides free and 0, 1, 2 or 3 intersections at the vertices; and 
those with no vertex intersections and one side intersection. 




The rest of the section describes how to divide a given 
triangle so that all the pieces are free triangles. We proceed 
step by step in terms of the number of free sides. 

Lemma 2.4. Every triangle with no free sides can be cut 
into seven free triangles. 

Proof. Each non-free side has one or two interior inter- 
sections with the conic. We draw the four cases and one 
solution for each (possible vertex intersections are irrelevant 
here, thus not drawn). All the small triangles can be proven 
free by noting that their sides are either free parts of the 
original sides, or segments connecting two intersections (thus 
free by part 2 of Remark 12. 2p . 




D 



The next step is to consider non-free triangles with one 
free side. We introduce another useful term. 

Definition 2.5. A triangle is almost free if exactly one 
of its sides is not free, and that side has only one intersection 
in its interior. 

Remark 2.6. There are four types of almost-free trian- 
gles, depending on the vertex intersections. 




Note that if a triangle is almost free and has no vertex in- 
tersections with the conic, then it is free (first case). 

Lemma 2.7. Every triangle with exactly one free side can 
be cut into five free or almost-free triangles, with zero or two 
of them being almost free. 



Proof. There are three cases depending on the number 
of interior intersections with the non-free sides: 2 + 2, 2 + 1 
and 1 + 1. The diagrams show how to cut the triangle in the 
three cases (possible vertex intersections, marked in white, 
make no difference). The numbers of free sides in each piece 
are indicated. As before, if a segment intersects a conic in 
its endpoints then it is free by part 2 of Remark 12.21 




In all cases, the only segment that may not be free is the 
lowest new one (dotted), and it can only have one interior 
intersection (by part 1 of Remark 12.21 since one of its end- 
points is in the conic already). D 

Next, we consider the case when two sides are free. 

Lemma 2.8. Every triangle with exactly two free sides can 
be cut into four free or almost-free triangles. At least one of 
them is free, except tf the original triangle is almost free. 

Proof. The non-free side has one or two interior inter- 
sections; in the former case, the triangle is almost free and 
we are finished. If it has two interior intersections, we join 
them with the opposite vertex and create two interior sides. 
There are three possibilities: 

1. If there are no interior intersections in the new sides, 
the three pieces are free. 

2. If one of the new sides has one interior point, we obtain 
one free triangle and two almost-free triangles. 

3. If both new sides have one interior point each, with 
one extra cut we obtain one free triangle and three 
almost-free triangles. 

The dashed lines indicate the partitions described above. 




Note that in the last two cases, cutting along the dotted line 
reduces the number of almost-free triangles by one, but it 
increases the total number of triangles. This might make a 
small difference in performance. D 

Lemma 2.9. Every almost-free triangle can be cut into 
four free triangles. 

Proof. If no vertex is in the conic, the triangle is already 
free (first case of Remark 12. 6p . If the vertex opposite to the 
non-free side belongs to the conic (third and fourth cases of 
Remark 1 2. 6 1) then the segment between them is free, and the 
triangle is cut into two free pieces. 

There remains only one case (see figure below): the conic 
intersects the triangle at two points, a vertex A and an in- 
terior point D of the side AB. The conic cannot intersect 
AB tangentially (otherwise it would have multiplicity inter- 
section > 3 with that line). Therefore it must enter the 
triangle through D and it can only exit through A. 



When CD is free, this segment cuts the triangle in two 
free pieces. However this is not true in general. 

Choose any point P in the conic and inside the triangle, 
with the property that BP and CP are free; then the original 
triangle is cut in four free pieces. 

C 




It suffices that the tangent to the conic at P leaves B and 
C on the same half-plane: the branch of the conic must 
then be contained in the other half-plane, thus BP and CP 
will be free. We offer three such points which are efficiently 
computable: the point whose tangent is parallel to BC; and 
the points at which the tangents pass through B or C. D 

Combining all the previous lemmas and counting the num- 
ber of pieces at each step, we obtain the following result. 

Proposition 2.10. Every triangle can be cut into eleven 
free triangles. 

Remark 2.11. It is possible to reduce the final number of 
free pieces to nine but one needs to use more often the re- 
course of finding tangency points in the come as in Lemma [279 
we chose the simpler approach. On the other hand, those 
points can be computed efficiently, which may make it attrac- 
tive to minimize the number of triangles in practice. Still, 
the integration time in each piece depends on the particular 
intersections. 

2.1 Degenerate conies 

We analyze now how to calculate the integral when / = 
is a degenerate conic. If it is empty, one point, or a double 
line, the integral is zero or the value on the full triangle. 

2.1.1 Two parallel lines 

If / = is two parallel lines, it can be converted by an 
affine transformation into x{x — 1) = 0. We can determine 
in which of the regions a:<0, 0<a:<l,l<x lies the 
image of each vertex by looking at their x-coordinates. 

1. If all three vertices are in one of the three regions, the 
integral is either the full triangle integral or zero; we 
can determine the sign of / in the triangle and use ^ . 

2. Otherwise, the triangle is split into two or three pieces 
(not necessarily triangular). The figure below depicts 
the possible cases. Once we have determined on which 
region(s) we must integrate (the middle strip or its 
complement), this can be done solely by adding and 
subtracting integrals of triangular pieces, which can be 
calculated using ([2]). 






2.7.2 Two crossing lines 

The conic can be transformed to the pair of hnes xy — 
0, what aUows us to quickly determine in which quadrants 
the vertices lie. The region on which to integrate is the 
intersection of the triangle and two opposing quadrants. 

1. If all vertices are in the same quadrant, the integral is 
the full triangle or zero. 

2. If all vertices are in two adjacent quadrants, the trian- 
gle is divided in two pieces, one of which is a triangle (or 
both, if a vertex lies in the limiting line) . The integral 
is that on the triangular piece, or the complementary. 

3. If the triangle is divided in three pieces by the conic, 
either all vertices are in different regions, or they are 
in two opposing regions. In any case, we can compute 
the integral on the relevant region by adding and sub- 
tracting integrals on triangles. 




Finally, if the triangle is divided in four pieces by the 
conic, there are two possible arrangements as well. Again, 
we can compute the integral by adding and subtracting 
triangles. 




For example, in the left figure, the union of the top- 
right and bottom-left regions of the triangle is a -I- 7 = 
{a + P + ^ + 5) — {l3 + ^) — (5 + ^) + 27, where the terms 
in parenthesis, as well as 7, are triangles. 

We can decide if we are in situation 1 or 2 by inspecting the 
signs of the coordinates of the transformed vertices. In order 
to differentiate situations 3 and 4 we use that the triangle is 
divided in four pieces if and only if the intersection of the two 
lines lies in its interior. This can be detected by calculating 
its barycentric coordinates as in the following algorithm. 

Algorithm 2.12. Determine if a point P is inside a tri- 
angle ABC. 

1. In the expression AP — aAB + {3 AC calculate a and 



det(A^,IC?) ^ det(Ii^,A^) 
det(AB,AC) 
where Aet{u,v) — uiV2 — U2«i. 



det 



(A§,IC?) 



2. If a, P > and a + p < 1 then P is contained in the 
triangle. 



3. If a ^ and f] e [0, 1]; or P ^ and a € [0, 1]; or 
a + /3 = 1 and a G [0, 1], then P is in the border of the 
triangle. 

4. Otherwise P is outside the triangle. 

In any case, for the final integral it is enough to add and 
substract several instances of Q, with no subdivisions other 
than the given by the lines of the conic. 

3. BASE CASE INTEGRATION 

In this section we describe how to detect the relative po- 
sition of the nondegenerate conic / = and a free triangle 
ABC, and compute the integral, in the five possible cases of 
free triangles. 

3.1 No intersections 

There are three possibilities: 



1- r C {/ > 0}: 
pleO 



the integral was computed in Exam- 



2. T C {/ < 0}: the integral is zero. 

3. / = is an ellipse contained in T. 

By inspecting / we can decide immediately whether f — 
is not an ellipse, from which we would deduce that we are in 
the first or second case. Then one can discern by evaluating 
the sign of / at some interior point of the triangle. 

On the other hand, if / = is an ellipse, we have to 
determine if any of the two shapes is contained in the other. 
We can do this by mapping the ellipse to the unit circle. 

Algorithm 3.1. Determine the relative position of the 
ellipse f — and the triangle ABC, and the correct domain 
of integration. 



1. Calculate an affine transformation 0: R — > 
sends / = mto x^ + y'^ = 1. Let P = (0, 0). 



M.-' that 



2. If d{(f>{A), P) < 1 then ABC is contained m the ellipse; 
evaluate the sign of f at some interior point of ABC 
to decide if the integral is the full triangle or zero. 

3. Otherwise, decide if P is mside the triangle A'B'C' := 
4>{ABC) with Alaonthm [2lM 

4. If P is in the triangle, then the ellipse is contained m 
it; evaluate the sign of f at (f>~^{P) to decide on which 
region to integrate. 

5. Otherwise, none of the shapes contains the other; eval- 
uate the sign of f at (fi~^{P) to decide if the integral is 
the full triangle or zero. 

The remaining computation is the integral of g when the 
ellipse / = is contained in the triangle. We show how to 
obtain a closed formula when {/ > 0} is the bounded region 
inside the ellipse; in the other case, the required integral is 
the difference of the full triangle integral and the former. 

Let (/p = 0"^ : R^ ^ R^ which sends the circle x^ -h y^ = 1 
to /. Then 



{/>()} 



g dxdy 



9{'P)\J{'P)\dxdy 



where D is the unit disc. Since ip is affine, |J(v3)| G R and 
5 := <?('P) is again a polynomial. Now, using polar coordi- 
nates, this is equal to 



g{r cos 6, r sin 9) ■ r dr\ d6 



\jm 



which is reduced to a linear combination of integrals of type 
J^'" cos'esm^ Ode. 

Alternatively, by Green's theorem the integral inside the 
ellipse is 



g{x,y)dxdy= G{x,y)dy 

J BE 



where 4^ = q. 

ox -^ 



3.2 One side intersection, no vertex intersec- 
tions 

This case is entirely similar to the previous one. 

3.3 One vertex intersection 

This case is even simpler: T is contained in {/ > 0} or 
{/ < 0}, we evaluate the sign of / at some interior point of 
the triangle in order to decide, and the integral will be that 
on the full triangle or zero. 

3.4 Two vertex intersections 

This case is more interesting. Either T is contained in 
one of the regions {/ > 0}, {/ < 0}, or it is divided in 
two regions by the conic. This can be discerned in the fol- 
lowing way: determine a segment which cuts the triangle in 
two (not necessarily triangular) pieces, separating the two 
relevant vertices, and count the number of intersections of 
that segment and the conic. Examples: the median of the 
side determined by the two vertices, or a suitable vertical or 
horizontal segment. If there are intersections, we are in the 
latter situation, otherwise evaluate the sign of / inside the 
triangle to separate the first two possibilities. 




Alternatively, convert the conic to a standard conic and 
check where the points lie after the transformation (see de- 
tails in the next subsection). 

If the triangle is divided in two regions by the curve, we 
really have to compute the integral on a region bounded 
by a conic arc and one or two segments. As usual we can 
consider only the former (the bottom region in the above 
picture), without loss of generality. How to determine the 
actual region of integration? The sign of / in the bottom 
region is the same as the sign of / in the middle point of the 
bottom side, for example. 

We can calculate the integral by afRnely transforming the 
conic into a standard conic: the circle x"^ -\- y^ — 1, the 
parabola y — x^ or the hyperbola xy = 1. 

1. Circle: the integral on the circular segment can be effi- 
ciently calculated as the integral on the circular sector 



minus the integral on the triangle determined by the 
segment and the center of the circle. 

Parabola: the integral after the transform is that on 
the region {y G [Iab{x),x'^],x £ [ai,6i]} where (ai, 02) 
and (61,62) are the images of the two intersection ver- 
tices, with ai < 61, and Iab{x) is the equation of the 
line through them. 




3. Hyperbola: similarly to the previous case, the integral 

can be calculated as that on the region {y £ [Iab{x), 1/x\,x £ 
[ai,6i]} if ai < 61 < 0, or {y e [l/x,lABix)],x G 
[ai,6i]} if < ai < 61. 

3.5 Three vertex intersections 

As in the previous case, either T is contained in one of 
the regions {/ > 0}, {/ < 0}, or it is divided in two regions 
by the conic. This time we use a different method to difi'er- 
entiate the three possibilities, since in the third one we also 
need to know which are the two vertices through which the 
conic enters the triangle. 

Since ellipses and parabolas define a convex region, a tri- 
angle with three vertices on such a curve cannot be divided 
by it. Thus, if the curve is of one of those types, it suf- 
fices once more to evaluate the sign of / in the triangle, and 
calculate the full triangle integral or return zero. 

If / = is a hyperbola, transform / = into xy = 1. This 
curve defines two convex regions, limited by the branches 
xy = l,x < and xy — l,x > 0. By inspecting the signs 
of the a;-coordinates of the (transformed) vertices, we can 
determine in which branch they are. 

1. If the three vertices are on the same branch of the 
hyperbola, the triangle is contained in {/ > 0} or {/ < 
0}, just determine the sign of / inside. 

2. Otherwise, two vertices lie on one branch and the third 
vertex lies on the other branch. The integral is calcu- 
lated as at the end of Section 13.41 




Note that the approach used in this case, namely the con- 
version to a standard conic in order to locate the vertices 
in relation to the curve, would have worked as well in Sec- 
tion l3.4l when we wanted to decide if the conic separates the 
triangle in two regions. This would amount to: 

1. Ellipse: convert to x^ + y^ = 1 and decide if the third 
vertex is inside or outside the unit circle. 



2. Parabola; convert to y = x^ and decide if the third 
vertex is above or below the parabola. 

3. Hyperbola: convert to xy = 1. If the two intersection 
vertices have different signs in their s-coordinates, the 
curve cannot separate the triangle. Otherwise, decide 
if the third vertex is in the convex region limited by 
the branch where the other two vertices are. 

4. THE ALGORITHM 

Algorithm l4.1l (next page) is a compilation of the steps de- 
scribed in the previous sections, so as to present an overview 
of the complete algorithm. Some case-by-case methods have 
not been explicitly written for brevity reasons. 

4.1 Practical considerations 

In relation to our implementation of this algorithm in 
MATLAB (almost complete as of May 2010) we would like 
to comment on numerical aspects that are not considered 
in our discussion above. First, several transformations sug- 
gested (Example 11.21 and the various transformations into 
standard conies from Section [3]) are a source of rounding 
errors because for small regions the scaling needed is very 
large. This problem can be solved by avoiding all scalings, 
i.e. restricting the transformations to rotations and transla- 
tions, not to a particular standard conic but to a member of 
some family of them. The result is a slight complication in 
the integration formulas, but nothing of concern in terms of 
efficiency. 

An additional problem is that in some cases (the calcu- 
lation suggested in Section 13.41 for the ellipse; Section 12. 1|) 
the sought integral is calculated as the difference of two easy 
integrals which may be orders of magnitude larger than the 
target, requiring much more precision in order not to lose 
significant digits. 

5. APPLICATION: AN OPTIMAL CONTROL 
PROBLEM 

Many technical processes are described by partial differ- 
ential equations. Here, it is important to optimize these 
processes. This leads to optimization problems in an infinite- 
dimensional setting. As an prototype, we consider the min- 
imization of a convex and quadratic functional subject to a 
linear elliptic partial differential equation and inequality con- 
straints on the control. Let us briefly introduce the optimal 
control problem we have in mind. 

Let n C R^ be a bounded domain with C'^-boundary F. 
For brevity, we will use ^ — {x,y) to denote points in R^. 
Let us introduce the following elliptic equation 



over all f £ L (0,) subject to the elliptic equation (|3]) and 
the control constraints 



-V ■ (DiOVuiO) + ciOuiO = Xn'fiO in Q, 
u{^) =0 on F. 



(3) 



Here, the control is denoted by /, while the solution u of 
this system is the corresponding state. Thanks to the as- 
sumptions below, for each control / £ L^{Q) there exists 
a unique response u £ i/o(fi), which is a weak solution of 
equation (|3}, see e.g. |2) Sect. 5.8]. The control acts on a 
compact polygonal subset fi' C fi. Now, we consider the 
control problem of minimizing 



fa < fiO < fb a.e. on Q. 



(5) 



That means, we want flnd a control / whose response u 
minimizes the distance to some desired state Ud. Let us 
denote this optimal control problem ©-© by (P). The set 
of admissible controls for (P) is given by 

Fad = {/ G L^{n) ■■ fa<f<fb a.c. on Q}. 

5.1 Existence and regularity of solutions 

Concerning the data of the state equation ((3)l, we make 
the following smoothness assumption on the data. 

Assumption 5.1. The coeffictents m the differential op- 
erator satisfy D £ C^'^{Q,) and c £ C^'^{Q,). Moreover, we 
assume that D{x) > Do > and c{x) > for all x £ Q. 

In order to obtain existence of solutions to (P) as well as 
a-priori discretization error estimates, we take the following 
assumptions on the data of the optimization problem. 

Assumption 5.2. We have a > 0, Ud e H^{^), and fa, 
/i, £ R with fa < fb a.e. on fi. 

Due to convexity, the problem under consideration is uniquely 
solvable, with solution denoted by (u*,/*). Moreover, the 
solution can be characterized by the following necessary op- 
timality conditions. These conditions are also sufficient since 
the optimal control problem is convex, see e.g. [51 Ch. 2]. 

Theorem 5.3. Let /* be the solution of (P) with associ- 
ated state u* . Then there exists an adjoint state p* £ H^{Q.) 
such that the adjoint equation 



-V • (D(OVp*(C)) + c(Op*(0 = {n-udm m n, 
and the variational inequality 



(6) 



(ar(c)+p*(e))(/(e)-r(0)de>o v/£f,, (?) 



are satisfied. Moreover, the following pointwise representa- 
tion of the optimal control holds 



r(C) = n/.,/.l(-V(C)) a.e. on n'. 



(8) 



Here, 'P[fa,fb](f) denotes the projection o/ / £ R on the in- 
terval [fa,fb]- 

Using the projection representation of the optimal control, 
we can conclude higher regularity of the solution: 



Theorem 5.4. Under the smoothness assumvtions \5.1\ and 
2J1 It holds u,p* £ H^{Q.), f £ H\Q,). 



J{f,u) 



{u(0-ud{,0fdi + 



.f\Od^ (4) 



Proof. Since we have p* £ H^{Q) by the previous theo- 
rem, the projection representation ^ implies that the opti- 
mal control has the same regularity /* £ H^{il). Then the 
right-hand sides of ^ and © are functions in H^{Q,). Stan- 
dard regularity results for elliptic partial differential equa- 
tions, e.g. [3 Thm. 8.13], yield u*,p* e H'^iQ). D 



Algorithm 4.1. Integrate a polynomial g{x,y) of degree 4 on the intersection of a triangle T and the region {/ > 0} 
determined by a quadratic polynomial f[x,y). 

1. If C :— {/ — 0} IS a degenerate conic, go to step 9. 

2. Calculate the intersections of C with each side of T . 

3. If all sides ofT are not free, let L :— {Ti, . . . , r„} be a list of free triangular pieces as in Lemma \2.4\ and go to step 6. 
4- Otherwise, use Lemma \2.7\ or Lemma \2.8\ to obtain a list L := {Ti, . . . , T„} of free or almost-free triangular pieces. 

5. For each triangle in L, if it is not free, substitute it m the list by the free pieces provided by Lemma \2.S\ 

6. Determine the type of C . 

7. Initialize S — 0. For each triangle Ti in L: 

7.1. Let Zi be the intersection of the border of Ti and C. 

1.2. If Zi — ^ or one non-vertex point: 

A. If C is an ellipse, use Alaorithm \3.1\ to know the relative position of C andT. 

I. If C is contained in T, determine the sign of f inside the ellipse. Let I be the integral of g on the bounded 

region inside C , or its complementary with respect to the full triangle, as needed, 
a. In any other case, determine the sign of f inside T. If it is positive, let I — Jj^, g, otherwise let I — 0. 

B. If C is not an ellipse, determine the sign of f inside T. If it is positive, let I — J J g, otherwise let 1 = 0. 

C. Add I to S. 

7.3. If Zi is one vertex: determine the sign of f in T. If it is positive let I = jjrp g, otherwise let 1 = 0. Add I to S. 

7.4. If Zi is two vertices: 

A. Calculate the number of intersections of C with the segment from the middle point of the two vertices to the 
third vertex. 

B. If there are none, determine the sign of f mside T. If positive, let I = JJ^, g, otherwise let 1 = 0. Add I to S. 

C. If there is one, determine which of the two regions is the correct one, by evaluating f in a suitable point. 

i. If C is an ellipse, transform it into x^ -\- y^ = 1. Calculate the integral on the circular segment. Let I be 

equal to that value or its complementary with respect to the full triangle, 
ii. If C is a parabola, transform it into y = x^ . Calculate the integral between the segment and the arc of 

parabola (the segment is always above). Let I be equal to that value or its complementary with respect to the 

full triangle. 
Hi. If C is a hyperbola, transform it into xy = 1. Calculate the integral between the segment and the arc of 

hyperbola (which one is above depends on which branch the vertices are in). Let I be equal to that value or 

its complementary with respect to the full triangle, 
iv. Add I to S. 

7. 5. If Zi is three vertices: 

A. If C is an ellipse or a parabola, determine the sign of f inside T. If it is positive, let I = J J g, otherwise let 
1 = 0. 

B. If C is a hyperbola, transform it into xy = 1 and determine in which branch does each vertex lie. 

i. All in one branch: determine the sign of f inside T. If it is positive, let I = JJj, g, otherwise let 1 = 0. 

ii. Two vertices A, B in one branch and the third vertex in the other branch: calculate the integral between 
the segment AB and the arc of hyperbola (which one is above depends on which branch the vertices are m). 
Determine the sign of f in the middle point of AB. If positive, let I be equal to the calculated integral; if 
negative, to its complementary with respect to the full triangle. 
C Add I to S. 

8. Output S and stop. 

9. Determine the type of degenerate conic. 

9.1. If C is empty, one point, or a double line, determine the general sign of f. If it is positive, let S = Jjrpg, otherwise 
let S = 0. Output S and stop. 

9.2. Otherwise, if C is two parallel lines, convert it to x — x = 0; if C is two crossing lines, convert it to xy = 0. 

9. 3. Determine the position of the vertices with respect to the lines by examining the coordinates of their images by the 
transformation. 

A. If all three vertices are in one of the regions, determine the sign of f inside T. If it is positive, let S = JJj,g, 
otherwise let S = 0. Output S and stop. 

B. Otherwise, determine the region(s) of integration by evaluating the sign of f at some vertex not on the conic. 
Write the region of integration as a sum of triangles with ±1 coefficients. Calculate the integral according to 
this. Output the result and stop. (A case by case method can be easily written.) 



5.2 Discretization and error estimate 

Now, we turn to the discretization of (P). To that end, 
let us introduce a family of quasi-uniform triangulations of 
Q, denoted by {Th}h>o- Each triangulation is assumed to 
exactly fit the boundary of O, such that Q = UreThT- This 
implies that elements of Tk lying on the boundary are curved. 
We further assume that for each T £Th there is a mapping 
"I>T mapping the standard simplex T to T. Moreover, we 
require that the intersection of every triangle T £ Th with 
the boundary of the control domain fl' is empty. That is, 
the boundary of Q' in fl is completely resolved by edges of 
triangles. 

With a triangulation we associate the following space of 
functions 

Vh = {ve C{n) : $t{v\t) G P2{f) VT G Th}, 

which implies that functions Vh G Vu are polynomials of 
degree 2 on each triangular element. Since fl' is a compact 
subset of Q, there is a mesh size /lo > such that all elements 
T £ Th with f2' n T 7^ are triangular. Hence, the above 
developed integration procedure can be applied for functions 
Vh G Vh with support in SI'. 

Then the discrete optimal control problem can be written 
as: minimize J{uh,fh) subject to Uh G 14, fh G Fad 

{DVuhVvh + cuhVh) d^ ^ fhVhd^ Vuh G Vh- 

n J Jiv 

(9) 
Note that we did not explicitly require fh to be in a finite- 
dimensional subspace. Nevertheless, if (ul^, fh) is a solution 
of the discrete problem, there exists a discrete adjoint state 
p*h G Vh satisfying 



//„ 



{DVplVvh + cplvh) di= I I (u'h~Ud)vhdS, Vuh G Vh 
n J Jn 

(10) 



and 



fh ^'Plfa,h] ("q^M 



(11) 



Due to this projection representation, the control is implic- 
itly discretized as the truncation of a function from the finite- 
dimensional space Vh- 

Theorem 5.5. Let (uh, fhjPh) be the solution of the dis- 
cretized optimality system ((9]) - (lllf) . Then there is a constant 
c > independent of the mesh size h such that 

Wfh - .f*llL2(n) + ll^h - ""'llfficn) + \\Ph -p*||H-i(n) < c K^ . 

Proof. Due to the approximation results of [I] Ch. 5.4], 
we have that the Assumption 2.4 in [3] is satisfied with Z = 
H''^{Q,) n HI{Q,) and convergence order h^ - Then the claim 
follows by a direct application of [31 Thm. 2.4]. D 

Known estimates for piecewise linear elements yield a conver- 
gence order of h^ only, compare [3]. In the two-dimensional 
case, i.e. fi C R'^, the number of unknowns A'^ — 2dimVii 
in the discretized problem is proportional to h~^ - Hence, 
our result implies that the approximation error is propor- 
tional to N~^^'^, whereas the use of linear polynomials only 
reduces the error like N~^ - This clearly shows that for opti- 
mal control problems as considered here, the use of piecewise 
quadratic approximations is preferable. 



5.3 Solution method 

In order to substitute fh in (fTCT)) by the projection ([TTJI. 
integrals 

faVhd^, // PhVhdS, 

have to be evaluated for piecewise quadratic polynomials 
Vh G Vh. This means, any solution method for the discretized 
problem encounters the difficulties of integrating over regions 
bounded by triangles and conies. 

The system consisting of the equation ((9)l- (|lip can be 
solved by means of a semi-smooth Newton method, see e.g. 
[3l. Within each step of the method, the non-smooth equa- 
tion (|lip is replaced by a linearized version 



X{-c«- 



ip!;e[/a,/b]} yjf' ^ ~P'' 



h 







on Q. , 



where p^ is the adjoint state given by the previous step, and 
Xa denotes the characteristic function of a set A- Multiply- 
ing this equation by a test function Vh G Vh and integrating 
on n', we obtain 

= // (fh - -P,^| Vhd£, 

iin'n{-Q-ip';e[/a,/i,]} V '^ / 



E 



fh Ph I Vhd^ 



for all Vh G Vh- Here, it is important to be able to evaluate 
the integrals 



IL 



Tn{-c-ip^elf^,ft]} 



'PlfaM ("qP'O "'^"^^ 



and 



// 

J Jrni- 



iphe[/a,/i,l} 



fhVhd^, 



which can be transformed to the type in the previous sec- 
tions. 

6. ACKNOWLEDGMENTS 

The authors would like to thank J. Schicho for his sug- 
gestions and the participants of the Rastenfeld workshop for 
their feedback. 

7. REFERENCES 

[1] S. C. Brenner and L. R. Scott. The mathematical theory 
of finite element methods, volume 15 of Texts in Applied 
Mathematics- Springer, New York, third edition, 2008. 

[2] D. Gilbarg and N. S. Trudinger. Elliptic partial 

differential equations of second order- Springer- Verlag, 
Berlin, 1983. 

[3] M. Hinze. A variational discretization concept in 

control constrained optimization: the linear-quadratic 
case. J- Computational Optimization and Applications, 
30:45-63, 2005. 

[4] F. Troltzsch. Optimale Steuerung partieller 

Differentialgleichungen- Vieweg, Wiesbaden, 2005. 

[5] R. J. Walker. Algebraic curves- Springer- Verlag, New 
York, 1978. Reprint of the 1950 edition. 



